home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 26
/
Cream of the Crop 26.iso
/
os2
/
octa209s.zip
/
octave-2.09
/
libcruft
/
odepack
/
prepj.f
< prev
next >
Wrap
Text File
|
1995-12-14
|
7KB
|
178 lines
SUBROUTINE PREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM,
1 F, JAC, IERR)
CLLL. OPTIMIZE
EXTERNAL F, JAC
INTEGER NEQ, NYH, IWM
INTEGER IOWND, IOWNS,
1 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
2 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
INTEGER I, I1, I2, IER, II, J, J1, JJ, LENP,
1 MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
DOUBLE PRECISION Y, YH, EWT, FTEM, SAVF, WM
DOUBLE PRECISION ROWNS,
1 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
DOUBLE PRECISION CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ,
1 VNORM
DIMENSION NEQ(*), Y(*), YH(NYH,*), EWT(*), FTEM(*), SAVF(*),
1 WM(*), IWM(*)
COMMON /LS0001/ ROWNS(209),
2 CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
3 IOWND(14), IOWNS(6),
4 ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH, MITER,
5 MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
C-----------------------------------------------------------------------
C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
C P = I - H*EL(1)*J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
C SUBJECTED TO LU DECOMPOSITION IN PREPARATION FOR LATER SOLUTION
C OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS DONE
C BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PREPJ USES THE FOLLOWING..
C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE).
C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE
C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
C OF P IF MITER IS 1, 2 , 4, OR 5.
C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
C EL0 = EL(1) (INPUT).
C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF
C P MATRIX FOUND TO BE SINGULAR.
C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
C (OR APPROXIMATION) IS NOW CURRENT.
C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
C MITER, N, NFE, AND NJE.
C-----------------------------------------------------------------------
NJE = NJE + 1
IERPJ = 0
JCUR = 1
HL0 = H*EL0
GO TO (100, 200, 300, 400, 500), MITER
C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
100 LENP = N*N
DO 110 I = 1,LENP
110 WM(I+2) = 0.0D0
CALL JAC (NEQ, TN, Y, 0, 0, WM(3), N)
CON = -HL0
DO 120 I = 1,LENP
120 WM(I+2) = WM(I+2)*CON
GO TO 240
C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
200 FAC = VNORM (N, SAVF, EWT)
R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC
IF (R0 .EQ. 0.0D0) R0 = 1.0D0
SRUR = WM(1)
J1 = 2
DO 230 J = 1,N
YJ = Y(J)
R = DMAX1(SRUR*DABS(YJ),R0/EWT(J))
Y(J) = Y(J) + R
FAC = -HL0/R
IERR = 0
CALL F (NEQ, TN, Y, FTEM, IERR)
IF (IERR .LT. 0) RETURN
DO 220 I = 1,N
220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
Y(J) = YJ
J1 = J1 + N
230 CONTINUE
NFE = NFE + N
C ADD IDENTITY MATRIX. -------------------------------------------------
240 J = 3
NP1 = N + 1
DO 250 I = 1,N
WM(J) = WM(J) + 1.0D0
250 J = J + NP1
C DO LU DECOMPOSITION ON P. --------------------------------------------
CALL DGEFA (WM(3), N, N, IWM(21), IER)
IF (IER .NE. 0) IERPJ = 1
RETURN
C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
300 WM(2) = HL0
R = EL0*0.1D0
DO 310 I = 1,N
310 Y(I) = Y(I) + R*(H*SAVF(I) - YH(I,2))
IERR = 0
CALL F (NEQ, TN, Y, WM(3), IERR)
IF (IERR .LT. 0) RETURN
NFE = NFE + 1
DO 320 I = 1,N
R0 = H*SAVF(I) - YH(I,2)
DI = 0.1D0*R0 - H*(WM(I+2) - SAVF(I))
WM(I+2) = 1.0D0
IF (DABS(R0) .LT. UROUND/EWT(I)) GO TO 320
IF (DABS(DI) .EQ. 0.0D0) GO TO 330
WM(I+2) = 0.1D0*R0/DI
320 CONTINUE
RETURN
330 IERPJ = 1
RETURN
C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
400 ML = IWM(1)
MU = IWM(2)
ML3 = ML + 3
MBAND = ML + MU + 1
MEBAND = MBAND + ML
LENP = MEBAND*N
DO 410 I = 1,LENP
410 WM(I+2) = 0.0D0
CALL JAC (NEQ, TN, Y, ML, MU, WM(ML3), MEBAND)
CON = -HL0
DO 420 I = 1,LENP
420 WM(I+2) = WM(I+2)*CON
GO TO 570
C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
500 ML = IWM(1)
MU = IWM(2)
MBAND = ML + MU + 1
MBA = MIN0(MBAND,N)
MEBAND = MBAND + ML
MEB1 = MEBAND - 1
SRUR = WM(1)
FAC = VNORM (N, SAVF, EWT)
R0 = 1000.0D0*DABS(H)*UROUND*DBLE(N)*FAC
IF (R0 .EQ. 0.0D0) R0 = 1.0D0
DO 560 J = 1,MBA
DO 530 I = J,N,MBAND
YI = Y(I)
R = DMAX1(SRUR*DABS(YI),R0/EWT(I))
530 Y(I) = Y(I) + R
IERR = 0
CALL F (NEQ, TN, Y, FTEM, IERR)
IF (IERR .LT. 0) RETURN
DO 550 JJ = J,N,MBAND
Y(JJ) = YH(JJ,1)
YJJ = Y(JJ)
R = DMAX1(SRUR*DABS(YJJ),R0/EWT(JJ))
FAC = -HL0/R
I1 = MAX0(JJ-MU,1)
I2 = MIN0(JJ+ML,N)
II = JJ*MEB1 - ML + 2
DO 540 I = I1,I2
540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC
550 CONTINUE
560 CONTINUE
NFE = NFE + MBA
C ADD IDENTITY MATRIX. -------------------------------------------------
570 II = MBAND + 2
DO 580 I = 1,N
WM(II) = WM(II) + 1.0D0
580 II = II + MEBAND
C DO LU DECOMPOSITION OF P. --------------------------------------------
CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
IF (IER .NE. 0) IERPJ = 1
RETURN
C----------------------- END OF SUBROUTINE PREPJ -----------------------
END